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LONG-TERM  GOALS 

The  overall  project  goal  has  been  to  test  the  feasibility  and  practicality  of  Bayesian  Hierarchical  Model 
(BHM)  methods  in  aspects  of  the  Mediterranean  Forecast  System  (MFS);  an  operational  ocean  data 
assimilation  and  forecast  system. 

OBJECTIVES 

Three  main  objectives  have  been  pursued  in  support  of  the  project  goal.  They  are: 

1.  a  surface  wind  BHM  (MFS-Wind-BHM)  to  drive  ensemble  ocean  data  assimilation  and  forecasts 
in  MFS; 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2Q  jq  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2010  to  00-00-2010 

4.  TITLE  AND  SUBTITLE 

Bayesian  Hierarchical  Models  to  Augment  the  Mediterranean  Forecast 
System 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Missouri, Department  of  Statistics, 146  Middlebush 

Hall, Columbia, MO, 6521 1 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

_ _ _  ABSTRACT 

18.  NUMBER  19a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  OS 

unclassified  unclassified  unclassified  Report  (SAR) 

7 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2.  a  time-  and  depth-dependent  background  error  covariance  BHM  (MFS-Error-BHM)  to  evolve  the 
background  error  covariance  matrix  in  13  sub-regions  of  the  MFS  forecast  domain;  and 

3.  a  BHM  to  demonstrate  super-ensemble  forecast  capabilities  (MFS-SuperEnsemble-BHM)  for 
ocean  applications. 

APPROACH 

MFS-Wind-BHM 

The  wind-driven  component  of  the  circulation  comprises  the  major  part  of  the  total  circulation  in  the 
Mediterranean  Sea  on  the  10-day  timescales  of  the  MFS  operational  forecasts.  MFS-Wind-BHM  is 
designed  to  perturb  the  surface  wind  forcing  for  data  assimilation  and  forecast  steps  of  the  MFS  ocean 
forecast  cycle  as  a  means  of  generating  an  ocean  ensemble  forecast.  The  surface  wind  perturbations  are 
constrained  to  be  within  the  uncertainties  prescribed  in  data  stage  and  process  model  distributions.  Data 
stage  distributions  are  constructed  from  measurement  error  models  for  the  QuikSCAT  surface  wind 
retrievals  and  from  analysis-forecast  error  statistics  for  the  ECMWF  surface  wind  fields.  The  process 
model  distribution  is  based  on  the  leading  order  terms  of  the  Rayleigh  Friction  Equations  (RFE;  Stevens 
et  al.,  2002).  The  purturbations  are  drawn  from  the  posterior  distribution  of  the  MFS-Wind-BHM  model 
for  the  14-day  sequential  data  assimilation  step  leading  to  ensemble  ocean  initial  conditions,  and  for  the 
10-day  forecast  step  generating  the  ensemble  forecast.  Ten  member  ensembles  are  used  in  the  pilot 
experiments.  Initial  condition  and  forecast  posterior  mean  fields  have  been  described  in  previous 
reports.  Initial  condition  and  forecast  uncertainty  are  quantified  by  the  spread  in  the  ensemble  fields  for 
the  ocean  state. 

MFS-Error-BHM 

The  approach,  recent  work  completed,  and  results  for  the  latest  implementation  of  MFS-Error-BHM  are 
presented  in  the  FY2010  annual  report  for  the  ONR  Physical  Oceanography  Program  BRC  project, 
“Bayesian  Hierarchical  Model  Characterization  of  Model  Error  in  Ocean  Data  Assimilation  and 
Forecasts”,  (Berliner,  Wikle,  Milliff;  co-PIs).  We  review  the  salient  points  in  abbreviated  sections  here 
and  refer  the  reader  to  that  report  for  more  detailed  discussion  and  figures. 

In  general,  the  dimension  of  the  full  state  {n  x  n)  error  covariance  matrix  Er  is  reduced  by  expanding  in 
a  set  of  / ix  p  basis  functions  (e.g.  EOFs,  where  n»  p)  as 

Er  =  <&Br<&' 

where  the  time  dependent  background  error  covariance  matrix  of  interest  is  Bt  of  dimension  p  x  p.  This 
matrix  is  not,  in  general,  diagonal  since  the  amplitudes  of  the  basis  functions  are  allowed  to  change  with 
time. 

We  define  an  error  process  model  as  et  —  Z/3r  +  r)t  where  Z  are  vertical  basis  functions  based  on 
multivariate  (i.e.  T  and  S)  profiles,  the  /3,  are  time-dependent  amplitudes  and  r)t  is  a  zero-mean 
Gaussian  error  term  to  account  for  additional  uncertainty  that  is  introduced  e.g.  in  the  dimension 
reduction. 

The  key  assumption  is  that  /3r  is  Gaussian,  with  zero  mean  and  variance  given  by  Bt  from  above.  We 
use  a  modified  Cholesky  decomposition  (Chen  and  Dunson,  2003)  to  write 

b,  =  A.r.r'rA, 
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where  At  is  a  diagonal  matrix  with  elements  proportional  to  the  standard  deviations  of  the  elements  of 
and  Tr  is  a  lower  triangular  matrix  associated  with  correlations  among  the  fit. 

At  the  next  level  of  the  Bayesian  hierarchy,  we  model  the  At  and  F,  as  autoregressive  time  series,  where 
the  autoregression  parameter  distributions  are  specified  at  still  lower  levels  of  the  BHM. 

Data  stage  distributions  are  built  from  misfit  dt  and  anomaly  qt  datasets  from  MFS  forecasts.  The  model 
misfits  are  forecast  differences  with  respect  to  in-situ  observations.  The  anomalies  are  departures  from 
the  model  “year  minus  day”  climatologies.  These  vectors  can  be  written 

dt  =  Ht(Xt\t_i)  —x°tbs  (1) 

dt  =  ~~  * 

where  Ht  is  the  operator  that  moves  the  forecast  Xt\t_i  to  the  observation  xf,s  locations  for  comparison, 
and  X,  is  the  climatology  value  for  the  model  state  variable  X. 

MFS-SuperEnsemble-BHM 

Super-ensembles  of  ocean  forecasts  in  the  Mediterranean  Sea  can  be  constructed  to  quantify  a 
multi-model  posterior  distribution  and  model  biases  for  target  ocean  processes.  Combining  ensembles 
from  several  models  is  a  challenging  problem  since  the  models  may  have  different  biases,  variabilities 
and  skills.  Adjusting  for  these  characteristics  and  managing  uncertainty  are  keys.  The  Bayesian  strategy 
suggested  in  Berliner  and  Kim  (2008)  offers  a  method  that  can  treat  these  issues.  In  this  approach 
ensembles  are  processed  much  like  observational  data  and  combined  with  a  prior  probability 
distribution  for  the  state  variables  of  interest. 

MFS-SuperEnsemble-BHM  is  our  implementation  of  the  Berliner  and  Kim  (2008)  methodology  for 
temperature  and  salinity  profile  evolution  (e.g.  T(z,t )  and  S(z-t  '))  at  two  locations  in  the  Rhodes  Gyre 
region  of  the  eastern  Mediterranean.  The  periods  of  interest  span  February-March  for  the  years 
2006-2009. 

We  have  assembled  ensembles  from  three  models:  NEMO,  OPA,  and  MedROMS.  Each  model  has  been 
forced  by  10  realizations  from  MFS-Wind-BHM,  as  well  as  by  the  ECMWF  analysis  winds  for  the 
periods  of  interest.  A  total  of  33  realizations  serve  as  our  data  stage  inputs.  The  MFS  analysis  fields, 
version  SYS3a2,  are  used  to  formulate  the  MFS-SuperEnsemble-BHM  prior  (or  process  model  stage 
distribution). 

Two  approaches  are  used  to  analyze  the  SYS3a2  output  to  form  a  prior: 

(1)  simple,  data  analytic  summaries  of  state  variable  means  and  covariances  serve  as  estimates  of 
prior  means  and  covariances;  and 

(2)  a  stochastic  time  series  model  (specifically,  a  multivariate  autoregressive  model  including 
parameters  for  two  temporal  lags)  is  fit  to  the  MedROMS  output  and  then  used  as  the  prior. 

WORK  COMPLETED 

MFS-Wind-BHM 

Two  manuscripts  (Milliff  et  al.,  2010;  and  Bonazzi  et  al.,  2010)  have  been  submitted,  reviewed,  and 
revised  for  publication  in  the  Quarterly  Journal  of  the  Royal  Meteorological  Society.  In  part  1  (Milliff 
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et  al.,  2010)  the  theory  and  implementation  of  MFS-Wind-BHM  is  described  in  detail.  An  appendix 
includes  the  term-by-term  definitions  of  the  probability  distributions  that  arise  at  each  level  of  the  model 
hierarchy.  In  part  2  (Bonazzi  et  al.,  2010),  the  impact  of  the  ensemble  forecast  methodology  based  on 
MFS-Wind-BHM  perturbations  is  documented.  Forecast  uncertainty  is  concentrated  at  the  ocean 
mesoscale;  it  is  shown  to  be  time  and  location  dependent. 

MFS-Error-BHM 

Wikle  has  been  computing  Bt  based  on  misfit  and  anomaly  data  stages  for  the  period  January-May 
2007,  for  the  MFS  sub-region  in  the  Gulf  of  Lions.  Wikle  supplies  vertical  EOFs  (Z)  and  Bt  to  Dr. 
Srdjan  Dobricic  at  MFS,  who  then  runs  twin  hindcast  experiments  with  the  operational  error  covariance 
and  the  Bt  from  MFS-Error-BHM  for  this  period.  Until  now,  there  have  been  “order-of-magnitude” 
differences  in  the  computed  error  covariances  between  the  two  methods. 

During  the  summer  visit  by  Prof.  Nadia  Pinardi  to  NWRA/CoRA  it  was  finally  determined  that  the 
regional  seasonal  cycle  removal  procedure  at  MFS  is  specific  and  should  be  implemented  in 
MFS-Error-BHM  as  well.  Preliminary  experiments  using  only  anomaly  (qt)  data  stage  inputs 
demonstrate  that  following  the  MFS  pre-processing  removes  the  discrepancy.  Figures  in  the  ONR 
annual  report  for  “Characterizing  Model  Error”  (see  Related  Projects)  demonstrate  that,  even  in  the 
absence  of  dt  data  stage  inputs,  the  forecast  impact  of  MFS-Error-BHM  is  neutral.  Experiments  are 
underway  now  to  introduce  dt  back  into  the  MFS-Error-BHM  and  quantify  forecast  impacts  at  MFS. 

MFS-SuperEnsemble-BHM 

We  have  assembled  all  needed  datasets  and  completed  algorithmic  development.  Archetypal  numerical 
codes  have  been  developed,  run,  and  assessed  for  effectiveness.  Currently,  those  codes  are  being  scaled 
up  to  the  multivariate  setting  of  our  S(z,t)  and  T{z,t)  datasets.  We  have  designed  experiments  to 
exemplify  and  assess  the  method,  including  the  role  of  the  prior  and  processing  the  Bayesian  output  for 
forecasting  purposes  as  well  as  model  assessment  (i.e.,  estimates  of  biases). 

A  major  paper  presenting  the  results  will  be  completed  and  submitted  in  December,  2010. 

While  the  code  development  for  MFS-SuperEnsemble  has  been  underway,  further  MedROMS 
development  and  validation  have  been  undertaken  by  Di  Lorenzo  and  colleagues  at  Georgia  Tech.  A 
MedROMS  ensemble  integration  (whole  basin,  3  members  so  far)  spanning  the  period  1963-2009  has 
been  produced  using  monthly  ECMWF  forcing. 

With  Dr.  Hazem  Nagy,  Di  Lorenzo  has  produced  an  early  draft  of  a  paper  on  the  climate  variability  of 
the  Mediterranean  eastern  basin  circulation.  The  ensemble  mean  reproduces  quite  well  the  satellite  data 
in  the  period  after  1990.  The  climate  simulations  also  capture  the  upper  circulation  inversion  that 
occurred  in  the  Ionian  sea  in  the  1990s.  Di  Lorenzo  et  al.  are  able  to  show  that  these  inversions  have 
occurred  in  the  past  on  multi-decadal  timescales.  These  mutli-decadal  fluctuations  associated  with  the 
upper  ocean  inversion  of  the  Ionian  circulation  are  captured  in  the  first  mode  of  sea-surface  height 
anomaly  variability  in  the  individual  realizations  and  in  the  ensemble  mean.  Di  Lorenzo  et  al  have  also 
isolated  an  atmospheric  bridge  from  the  Pacific  to  the  Mediterranean  Sea  in  their  analyses.  It  appears 
that  the  second  mode  of  variability  in  the  eastern  basin  circulation  integrates  the  effects  of  interannual 
tropical  Pacific  variability;  so  much  so  that  the  principal  component  of  the  eastern  basin  circulation 
mode  reproduces  indices  of  Pacific  climate  (i.e.  ENSO  and  PDO)  with  correlation  R  =  0.86.  The 
Georgia  Tech  group  is  exploring  mechanisms  and  implications  for  the  deeper  circulation  and  on  deep 
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water  formation  in  the  Mediterranean.  Interestingly,  the  effect  of  the  NAO  does  not  appear  to  dominate 
in  the  eastern  basin. 

While  the  MedROMS  climate-scale  research  is  not  a  task  directly  associated  with  any  of  the  MFS  BHM 
objectives,  it  does  demonstrate  model  stability  and  utility  in  this  newest  component  of  the 
MFS-SuperEnsemble-BHM  data  stage  inputs.  As  noted,  a  climate-scale  manuscript  is  in  preparation 
and  ONR  funding  will  be  acknowledged. 

Relevant  Presentations 

(Berliner,  Milliff,  Pinardi,  Wikle)  Informal  presentations  and  discussions  at  the  annual  “All-Hands” 
project  meeting  at  NWRA/CoRA,  August,  2010. 

(Milliff,  Wikle;  session  co-conveners)  Probabilistic  Models  in  Ocean  Sciences:  Applications  in  Data 
Assimilation,  Coupled  Ecosystem  Models  and  Air-Sea  Interaction  Studies,  American  Geophysical 
Union,  Ocean  Sciences  Meeting,  Portland,  OR,  February,  2010. 

(Berliner)  Combining  Models  and  Data:  The  Bayesian  Approach  to  Modeling  and  Prediction,  Invited 
Talk,  AGU  Ocean  Sciences  Meeting,  Portland,  OR,  February,  2010. 

(Pinardi)  A  new  method  for  ocean  ensemble  forecasting  with  quantification  of  wind  uncertainties. 
Invited  talk.  AGU  Ocean  Sciences  Meeting,  Portland,  OR,  February,  2010. 

(Milliff,  Pinardi,  Wikle,  Berliner,  Bonazzi)  Process  model  considerations  for  a  surface  wind  Bayesian 
hierarchical  model.  Poster,  AGU  Ocean  Sciences  Meeting,  Portland,  OR,  February  2010. 

(Milliff)  Estimating  semivariograms  to  build  covariance  matrices  for  J .  Workshop  on  the  ROMS  4D-Var 
Data  Assimilation  Systems  for  Advanced  ROMS  Users,  University  of  California,  Santa  Cruz,  July  2010. 

(Wikle)  A  hierarchical  approach  to  motivate  spatio-temporal  statistical  models.  Institute  for  Pure  and 
Applied  Mathematics  (IPAM),  UCFA,  Eos  Angeles,  CA,  May  25,  2010. 

(Wikle)  Bayesian  hierarchical  models  to  augment  the  Mediterranean  forecast  system.  Invited  talk.  Iowa 
State  University.  Ames,  IA,  October  15,  2009. 

(Wikle)  Don’t  forget  the  process!  Using  scientific  process  knowledge  to  motivate  spatio-temporal 
models.  Invited  talk.  SAMSI  Program  on  Space-Time  Analysis  for  Environmental  Mapping, 
Epidemiology  and  Climate  Change,  Opening  Workshop,  RTP,  North  Carolina,  September  14,  2009. 

(Wikle)  A  class  of  nonlinear  spatio-temporal  dynamic  models.  Invited  Talk,  Joint  Statistics  Meetings, 
Washington,  DC,  August  4,  2009. 

RESULTS 

MFS-Wind-BHM 

We  are  awaiting  word  from  QJRMS  on  the  status  of  the  revised  manuscripts;  Milliff  et  al.,  2010  (part  I) 
and  Bonazzi  et  al.,  2010  (part  II). 

MFS-Error-BHM 

Embedded  scales  in  the  error  covariance  estimations  of  ocean  forecast  systems  act  to  rescale  the  error 
covariance  magnitudes.  Anomaly  data  stage  inputs  are  probably  not  sufficient  to  represent  abrupt 
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regime  shifts  in  the  ocean  state  and  provide  added  value  over  operational  methods  using  a  fixed, 
seasonal,  background  error  covariance.  Experiments  adding  misfit  data  stage  inputs  will  be  useful  in 
modelling  error  covariance  response  to  ocean  regime  shifts  in  the  Mediterranean  on  sub-seasonal 
timescales. 

Under  seperate  ONR  funding,  the  MFS-Error-BHM  methodology  is  being  adapted  to  the  California 
Current  System  and  the  Regional  Ocean  Model  System,  four-dimensional  variational  data  assimilation 
and  forecast  tools  developed  by  Prof.  A.  Moore  and  colleagues. 

IMPACT/APPLICATIONS 

MFS-Wind-BHM 

BHM  methods  are  proving  feasible  and  practical  in  an  array  of  ocean  data  assimilation  and  forecast 
applications.  The  papers,  Milliff  et  al.  (2010)  and  Bonazzi  et  al.  (2010)  will  serve  as  detailed  generic 
references  for  this  methodology  in  ocean  and  atmospheric  sciences. 

MFS-Error-BHM 

Refining  estimates  of  the  time-dependent  changes  in  ocean  states  and  forecast  uncertainty  across  regime 
shifts  adds  value  to  ocean  forecast  system  output. 

MFS-SuperEnsemble-BHM 

First  steps  in  developing  super-ensemble  ocean  forecast  techniques  for  targeted  fields  and/or  processes 
will  provide  for  feasibility  and  practicality  assesments  of  this  popular,  but  poorly  understood, 
methodology.  Model  bias  distributions  will  demonstrate  a  new  means  of  comparing  forecast  tools  that 
are  being  developed  in  the  community  (i.e.  ROMS,  OPA,  NEMO). 

TRANSITIONS 

Informal  communications  continue  with  scientists  in  the  Ocean  Modelling  branch  of  the  Naval 
Research  Laboratory,  Bay  St.  Louis,  MI. 

RELATED  PROJECTS 

“Estimating  Ecosystem  Model  Uncertainties  in  Pan-Regional  Syntheses  and  Climate  Change  Impacts 
on  Coastal  Domains  of  the  North  Pacific  Ocean”,  NSF  US  Globec  Program,  October  2008  -  October 
2011. 

“Bayesian  Hierarchical  Model  Characterization  of  Model  Error  in  Ocean  Data  Assimilation  and 
Forecasts”,  ONR  Physical  Oceanography  Program;  Milliff,  Berliner,  Wikle,  co-PIs. 

“Quantifying  the  Amplitude,  Structure  and  Influence  of  Model  Error  during  Ocean  Analysis  and 
Forecast  Cycles”,  ONR  Physical  Oceanography  Program,  A.  Moore  (PI). 
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